#include <stdio.h>
#include <vector>
//#define CUDA_ENABLE
#include "cuda_lib.h"
#include "math.h"
#include "glshower.h"
#include "ShapeObject.h"
#pragma once

#define PLANET_RADIUS 6371e3 //地球半径6371千米
#define ATMOS_RADIUS 6471e3 //大气半径6471千米

// and the heights (how far to go up before the scattering has no effect)
#define HEIGHT_RAY 8e3 //rayleigh height（海拔8千米）
#define HEIGHT_MIE 1.2e3 //and mie（海拔1.2千米）
#define HEIGHT_ABSORPTION 30e3 //at what height the absorption is at it's maximum
#define ABSORPTION_FALLOFF 4e3 //how much the absorption decreases the further away it gets from the maximum height
#define PRIMARY_STEPS 12 
#define LIGHT_STEPS 4
constexpr float G = 0.7; //mie scattering direction, or how big the blob around the sun is
constexpr float GG = G * G; //mie scattering direction, or how big the blob around the sun is

// scattering coeffs
#define RAY_BETA vec3(5.5e-6, 13.0e-6, 22.4e-6) //rayleigh, affects the color of the sky。瑞利散射系数
#define MIE_BETA vec3(21e-6,21e-6,21e-6) //mie, affects the color of the blob around the sun。米散射系数
#define AMBIENT_BETA vec3(0.0) //ambient, affects the scattering color when there is no lighting from the sun
#define ABSORPTION_BETA vec3(2.04e-5, 4.97e-5, 1.95e-6) //what color gets absorbed by the atmosphere (Due to things like ozone)
#define G 0.7 //mie scattering direction, or how big the blob around the sun is

struct AtmosphereParamaters
{
	float view_pos[3];
	float light_dir[3];
};
__DEVICE__ void render_scene(const vec3& view_pos, const vec3& view_dir, const vec3& light_dir,vec3& scene_color, float& scene_depth)
{
	scene_color = vec3(0.0, 0.0, 0.0);
	scene_depth = 1e12;

	//画太阳
	scene_color = dot(view_dir, light_dir) > 0.9998 ? 3.0 : 0.0;

	uint8_t planet_data[sizeof(SphereShape)];
	auto planet = (SphereShape*)planet_data;
	planet->center = vec3(0,0,0);
	planet->radius = PLANET_RADIUS;
	IsectData planet_intersect;
	if (planet->intersect(view_pos, view_dir, planet_intersect))//if (0.0 < planet_intersect.t1)		
	{
		scene_depth = d_max(planet_intersect.t0, 0.0f);//交点中近处的那个点
		vec3 sample_pos = view_pos + (view_dir * planet_intersect.t0) - planet->center;
		vec3 surface_normal = normalize(sample_pos);
		//planet color
		scene_color = vec3(0.0, 0.25, 0.05);
		// get wether this point is shadowed, + how much light scatters towards the camera according to the lommel-seelinger law
		vec3 N = surface_normal;
		vec3 V = -view_dir;
		vec3 L = light_dir;
		float dot_nv = d_max<float>(1e-6, dot(N, V));
		float dot_nl = d_max<float>(1e-6, dot(N, L));
		float shadow = dot_nl / (dot_nl + dot_nv);

		scene_color *= shadow;
	}
}
/// <summary>
/// calculate the density of the particles (both for rayleigh and mie)大气密度值（也可叫作光衰减率,与高度有关）
/// and the absorption density. this is for ozone, which scales together with the rayleigh, 
/// but absorbs the most at a specific height, so use the sech function for a nice curve falloff for this height
/// clamp it to avoid it going out of bounds. This prevents weird black spheres on the night side
/// </summary>
/// <param name="altitude">海拔</param>
/// <returns>vec3(rayleigh_scater, mie_scater, absorb)</returns>
__DEVICE__ vec3 atmosphere_density(float altitude)
{
	const float denom = (HEIGHT_ABSORPTION - altitude) / ABSORPTION_FALLOFF;
	const float rayleigh_scater = exp(-altitude / HEIGHT_RAY);
	return vec3(rayleigh_scater, exp(-altitude / HEIGHT_MIE), (1.0 / (denom * denom + 1.0)) * rayleigh_scater);
}
__DEVICE__ float phase_rayleigh(float dot_view_light_sq)
{
	return 3.0 / 50.2654824574 * (1.0 + dot_view_light_sq);//50.2654824574:(16 * pi)相位函数;
}

__DEVICE__ float phase_mie(float dot_view_light,float dot_view_light_sq, float g, float gg)
{
	return 3.0 / 25.1327412287 * ((1.0 - gg) * (dot_view_light_sq + 1.0)) / (pow(1.0 + gg - 2.0 * dot_view_light * g, 1.5) * (2.0 + gg));
}

template<class T>
__DEVICE__ const T& d_read_param(const uint8_t* params, uint32_t& offset)
{
	auto ptr = (T*)(params + offset);
	offset += sizeof(T);
	return *ptr;
}

__DEVICE__ void atmosphere(color* image, const uint32_t width, const uint32_t height, uint32_t x, uint32_t y, uint32_t z);

#ifdef CUDA_ENABLE
__GLOBAL__ void atmosphere(color* image, const uint32_t width, const uint32_t height)
{
	uint32_t x, y, z;
	d_coord(x, y, z);
	atmosphere(image, width, height, x, y, z);
}
#endif

__CONSTANT__ AtmosphereParamaters atmosphere_paramaters;

__DEVICE__ void atmosphere(color* image, const uint32_t width, const uint32_t height, uint32_t x, uint32_t y, uint32_t z)
{
	if (x >= width || y >= height)
		return;

	vec3 fragColor(0, 0, 0);

	const vec3 PLANET_CENTER(0, 0, 0);
	const auto view_pos = *(const vec3*)atmosphere_paramaters.view_pos;// d_read_param<vec3>(params, params_top);// vec3(params[0], params[1], params[2]) - PLANET_CENTER;//相对位置
	const auto light_dir = *(const vec3*)atmosphere_paramaters.light_dir;// normalize(d_read_param<vec3>(params, params_top));

	const vec3 light_intensity = vec3(40.0);
	const vec3 view_dir = d_raydir(width, height, x, y, z);
	vec3 scene_color(0,0,0);
	float scene_depth{0};
	render_scene(view_pos, view_dir, light_dir, scene_color, scene_depth);
	fragColor = scene_color;
	//
	uint8_t atmosphere_data[sizeof(SphereShape)];
	auto atmosphere = (SphereShape*)atmosphere_data;
	atmosphere->center = vec3(0, 0, 0);
	atmosphere->radius = ATMOS_RADIUS;
	IsectData atmos_intersect;
	if (atmosphere->intersect(view_pos, view_dir, atmos_intersect)
		&& atmos_intersect.t0 < atmos_intersect.t1)
	{
		// prevent the mie glow from appearing if there's an object in front of the camera
		const bool allow_mie = scene_depth > atmos_intersect.t1;//镜头像素上没物体时
		const float mu = dot(view_dir, light_dir);
		const float mumu = mu * mu;
		const float phase_ray = phase_rayleigh(mumu);
		const float phase_m = allow_mie ? phase_mie(mu, mumu, G, GG) : 0.0;

		vec3 total_ray = vec3(0.0); // for rayleigh
		vec3 total_mie = vec3(0.0); // for mie

		// prevent the mie glow from appearing if there's an object in front of the camera
		atmos_intersect.t0 = d_max(atmos_intersect.t0, 0.f);
		atmos_intersect.t1 = d_min(atmos_intersect.t1, scene_depth);
		const float step_size_i = (atmos_intersect.t1 - atmos_intersect.t0) / float(PRIMARY_STEPS);//
		// initialize the optical depth. This is used to calculate how much air was in the ray
		vec3 opt_i = vec3(0.0);
		for (int i = 0; i < PRIMARY_STEPS; ++i)
		{
			const vec3 view_sample_pos = view_pos + view_dir * (i + 0.5f) * step_size_i;
			const float altitude_i = view_sample_pos.length() - PLANET_RADIUS;
			const vec3 dss = atmosphere_density(altitude_i);
			const vec3 density = atmosphere_density(altitude_i) * step_size_i;//密度值*分段距离
			//if (d_ispixel(width, height, 400, 200))
			//	printf("%f,%f,%f  %f  %f\n", dss.x, dss.y, dss.z, altitude_i, step_size_i);

			// Add these densities to the optical depth, so that we know how many particles are on this ray.
			opt_i += density;//视线上的空气密度积分

			//
			vec3 opt_light(0, 0, 0);
			IsectData light_atmo_intersect;
			atmosphere->intersect(view_sample_pos, light_dir, light_atmo_intersect);
			const float step_size_light = (light_atmo_intersect.t1 - light_atmo_intersect.t0) / LIGHT_STEPS;
			for (int light_step_idx = 0; light_step_idx < LIGHT_STEPS; ++light_step_idx)
			{
				vec3 light_sample_pos = view_sample_pos + light_dir * (light_step_idx + 0.5) * step_size_light;
				const float altitude_light = light_sample_pos.length() - PLANET_RADIUS;//海拔
				const vec3 density_light = atmosphere_density(altitude_light) * step_size_light;
				opt_light += density_light;
			}
			vec3 attn = exp(-RAY_BETA * (opt_i.x + opt_light.x) - MIE_BETA * (opt_i.y + opt_light.y) - ABSORPTION_BETA * (opt_i.z + opt_light.z));

			// accumulate the scattered light (how much will be scattered towards the camera)
			total_ray += attn* density.x;//density
			total_mie += attn* density.y;
		}

		//不透明度（exp(-(raylei scatter + mei scatter + absorb))
		const vec3 opacity =exp(-(RAY_BETA * opt_i.x + MIE_BETA * opt_i.y + ABSORPTION_BETA * opt_i.z));//视线方向上的不透明度
		
		vec3 in_scatter = (RAY_BETA * total_ray * phase_ray + MIE_BETA * total_mie * phase_m + AMBIENT_BETA * opt_i) * light_intensity;
		//vec3 in_scatter = (RAY_BETA * total_ray * phase_ray);
		//in_scatter += MIE_BETA * total_mie * phase_m;
		//in_scatter += AMBIENT_BETA * opt_i;
		//in_scatter *= light_intensity;
		
		fragColor = in_scatter + scene_color * opacity;// + 背景色 * 视线上的不透明度
	}//end-if
	fragColor = vec3(1.0) - exp(-fragColor);

	d_pixel(image, width, height, x, y, z) = {
		static_cast<uint8_t>(fragColor.x * 255),
		static_cast<uint8_t>(fragColor.y * 255),
		static_cast<uint8_t>(fragColor.z * 255),
		255 };
}